Generamos datos random de la distribución, poniendo atención al tercer parámetro, que es el valor de alfa
Sacamos los datos a graficar aplicando las operaciones pertinentes
Hacemos los datasets para graficar la relación lineal utilizando el valor de alfa como la pendiente
Por últmo, graficamos esta relación lineal
Primero se especifica la media y la mediana reportada por el U.S. Census Buerau. En este ejemplo se utiliza la del año 2008. La informaci?n puede encontrarse en https://www.census.gov/data/tables/time-series/demo/income-poverty/historical-income-households.html Las cantidades estan en miles de dolares.
Utilizando las relaciones entre la media y la varianza de la normal asociada calculamos los parametros para dicha normal
Restamos la media a la media mas la varianza para poder obtener la varianza solita y procedemos a sacar la raiz para obtener la desviacion estandar
Se hace un array con los cuantiles para los cuales se quiere estimar el ratio de desigualdad
quant = c(0, .1, .2, .3, .4, .5, .6, .7, .75, .8, .85, .9, .95, .98, .99, .995, .999, .9995, .9999, .99999)Utilizando la media y la desviaci?n est?ndar que ya obtuvimos, podemos entonces calcular los ingresos por household que corresponden a cada cuantil
Posteriormente, utilizando la relaci?n entre la distribuci?n de hogares por ingreso y la distribuci?n del ingreso total, podemos calcular los par?metros de la distribuci?n del ingreso total de la siguiente forma
Con estos parámetros, podemos estimar el porcentaje del ingreso total por ingreso
Se calcula entonces el área que hay a partir de cada punto
en quant_inverso se tiene el porcentaje de hogares que tienen un ingreso mayor a los puntos en ing_por_household y en int_total_per_inverso esta el porcentaje total del ingreso que acaparan los hogares que tienen un ingreso mayor a los puntos en ing_per_household
Para calcular el coeficiente relativo, unicamente dividimos el porcentaje del ingreso total acaparado por hogares con un ingreso mayor a un punto determinado entre el porcentaje de (el número) de casas que tienen un ingreso mayor al mismo punto.
Por ultimo se pone todo en un dataframe para que tenga el formato de tabla
tabla_final = data.frame("cuantil" = quant,
"Ingreso en K dolares" = ing_por_household,
"Ingreso Total hasta k dolares" = ing_total_per,
"P(X>x)" = quant_inverso,
"P(Y>y)" = int_total_per_inverso,
"ratio"= coef_rel)
print(tabla_final)## cuantil Ingreso.en.K.dolares Ingreso.Total.hasta.k.dolares P.X.x.
## 1 0.00000 0.00000 0.00000000 1.00000
## 2 0.10000 14.43286 0.01218988 0.90000
## 3 0.20000 22.11017 0.03505797 0.80000
## 4 0.30000 30.07204 0.06759562 0.70000
## 5 0.40000 39.11058 0.11068622 0.60000
## 6 0.50000 50.00000 0.16613799 0.50000
## 7 0.60000 63.92133 0.23693621 0.40000
## 8 0.70000 83.13370 0.32810958 0.30000
## 9 0.75000 96.15559 0.38397786 0.25000
## 10 0.80000 113.07013 0.44910674 0.20000
## 11 0.85000 136.57670 0.52666683 0.15000
## 12 0.90000 173.21579 0.62248424 0.10000
## 13 0.95000 246.35508 0.75026183 0.05000
## 14 0.98000 366.21264 0.86086397 0.02000
## 15 0.99000 476.99674 0.91257891 0.01000
## 16 0.99500 607.52408 0.94589485 0.00500
## 17 0.99900 1000.37074 0.98302616 0.00100
## 18 0.99950 1214.78132 0.98985623 0.00050
## 19 0.99990 1840.43491 0.99701548 0.00010
## 20 0.99999 3124.42020 0.99950851 0.00001
## P.Y.y. ratio
## 1 1.0000000000 1.000000
## 2 0.9878101217 1.097567
## 3 0.9649420294 1.206178
## 4 0.9324043786 1.332006
## 5 0.8893137759 1.482190
## 6 0.8338620125 1.667724
## 7 0.7630637877 1.907659
## 8 0.6718904199 2.239635
## 9 0.6160221442 2.464089
## 10 0.5508932573 2.754466
## 11 0.4733331743 3.155554
## 12 0.3775157585 3.775158
## 13 0.2497381734 4.994763
## 14 0.1391360276 6.956801
## 15 0.0874210871 8.742109
## 16 0.0541051484 10.821030
## 17 0.0169738402 16.973840
## 18 0.0101437723 20.287545
## 19 0.0029845247 29.845247
## 20 0.0004914939 49.149394
Simulamos primero datos aleatorios de una distribución Singh - Maddala, censurando artificialmente los datos que sean mayores a 3. De forma preeliminar, calculamos el porcentaje de los datos que estamos censurando y cual es la media real de los agentes con un ingreso mayor a 3
ingresos <- rsinmad(100000, 1.14, 1.75, 2.07)
top_code <- 3
prop_real <- sum(ingresos>top_code)/100000
real_mean <- mean(ingresos[ingresos>top_code])Procedemos a hacer la censura de los datos
ingresos_censurado <- rep(0,100000)
ingresos_censurado[ingresos < top_code] <- ingresos[ingresos < top_code]
ingresos_censurado[ingresos >= top_code] <- top_codeSuponemos entonces que la distribución de los ingresos para los datos arriba del cuantil 95% sigue una distribución Pareto Tipo I, y usando el estimador propuesto por Armour, Burkhauser y Larrimore estimamos el parámetro \(\alpha\).
Una vez que tenemos dicho parámetro, utilizamos la propiedad de la distribución pareto para calcular la media del ingreso de los agentes con un ingreso mayor a 3 y utilizamos esta media para imputar los datos que habían sido censurados. Esta nueva serie conservará mejor ciertas propiedades de los datos sin censura, como la media.
cutoff <- quantile(ingresos, .95)[[1]]
M = sum(ingresos_censurado>cutoff & ingresos_censurado<top_code)
Te = sum(ingresos_censurado==3)
alpha_est = M/(Te*log(top_code) + sum(log(ingresos_censurado[ingresos_censurado>=cutoff & ingresos_censurado<top_code])) - (M+Te)*log(cutoff))
invertido = alpha_est/(alpha_est-1)
estimated_mean <- top_code*invertido
## Media estimada
estimated_mean## [1] 4.562036
## [1] 4.28969
Notamos que la media estimada y la real son muy similares.
x <- rnorm(100000,0,1)
y <- rnorm(100000,0,1)
t <- x[x<1.5]
c <- y
c[y>=1.5] <- 1.5
x <- x+1
ggplot()+
geom_density(aes(x=t), fill = "#f04546", size = 1)+
geom_density(aes(x=c), fill = "#6db823", size = 1) +
geom_density(aes(x=x), fill = "#5b5c5a", size = 1, alpha=.4) +
geom_vline(xintercept = 1.5, linetype = 2, color = "#3591d1", size = 1) +
xlab("x") + ylab("Densidad")+
theme_bw()## ENIGH
## Homogeneizamos los nombres de las columnas de interes
concentradohogar2012 <- read_csv(here("data", "concentradohogar2012.csv"))
colnames(concentradohogar2012)[colnames(concentradohogar2012) == "factor_hog"] <- "factor"
concentradohogar2014 <- read_csv(here("data", "concentradohogar2014.csv"))
colnames(concentradohogar2014)[colnames(concentradohogar2014) == "factor_hog"] <- "factor"
concentradohogar2016 <- read_csv(here("data", "concentradohogar2016.csv"))
concentradohogar2018 <- read_csv(here("data", "concentradohogar2018.csv"))
## Datos anonimos del SAT
sat2012 <- read_delim(here("data", "RUIDO_FINAL_PF_2012.tab"), "\t", escape_double = FALSE, trim_ws = TRUE)Creamos histogramas con diferentes parámetros para familiarizarnos con los datos.
df <- tibble(ing_cor = c(concentradohogar2012$ing_cor,
concentradohogar2014$ing_cor,
concentradohogar2016$ing_cor,
concentradohogar2018$ing_cor),
factor = c(concentradohogar2012$factor,
concentradohogar2014$factor,
concentradohogar2016$factor,
concentradohogar2018$factor),
año = c(rep("2012", length(concentradohogar2012$ing_cor)),
rep("2014", length(concentradohogar2014$ing_cor)),
rep("2016", length(concentradohogar2016$ing_cor)),
rep("2018", length(concentradohogar2018$ing_cor)))
)
histograma <- ggplot(df, aes(x = ing_cor, weight = factor)) +
geom_histogram(bins = 60) +
facet_wrap(~año) +
scale_x_continuous(labels = scales::comma) +
scale_y_continuous(labels = scales::comma) +
xlab("Ingreso Corriente") + ylab("Frecuencia") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
histograma_zoomed <- ggplot(df, aes(x = ing_cor, weight = factor)) +
geom_histogram(bins = 600) +
facet_wrap(~año) +
scale_x_continuous(labels = scales::comma) +
scale_y_continuous(labels = scales::comma) +
coord_cartesian(xlim=c(0,3000000)) +
xlab("Ingreso Corriente") + ylab("Frecuencia") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
histograma_log <- ggplot(df, aes(x = ing_cor, weight = factor)) +
geom_histogram(bins = 30) +
facet_wrap(~año) +
scale_x_continuous(trans = "log", labels = scales::comma) +
scale_y_continuous(labels = scales::comma) +
xlab("Ingreso Corriente (log)") + ylab("Frecuencia") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))Ya que los histogramas no nos fueron de gran utilidad, nos pasamos a investigar qué pasa en cada decil de la distribución de los ingresos. Además, buscamos ver si existe consistencia entre cada año de aplicación de la ENIGH.
ingreso_medio_percentil <- function(datos, n_percentil) {
# datos son las base de datos del censo, n_percentil es el numero de divisiones (e.g.: 10 retorna deciles, 100 retorna percentiles)
## Datos sorteados por ingresos trimestales corrientes más indice de los hogares representados observados
sorted_ing <- datos %>% select(ing_cor, factor) %>%
arrange(ing_cor) %>%
mutate(cumulative_obs = cumsum(factor))
## Total de hogares representados (ver definicion de factor en documentacion de la ENIGH)
total_hogares <- sum(datos$factor)
## Tamaño por decil truncado para no incluir decimales
tam_decile <- trunc(total_hogares/n_percentil)
## Primeros n-1 deciles
lim_inf <- 1
lim_sup <- tam_decile
ing_med_de <- c()
for (i in 1:(n_percentil-1)) {
ing_med_de[i] <- sorted_ing %>% filter(cumulative_obs >= lim_inf, cumulative_obs <= lim_sup) %>%
summarise(mean = sum(ing_cor*factor)/tam_decile) %>%
pull()
## Limites del decil
lim_inf <- lim_inf + tam_decile
lim_sup <- lim_sup + tam_decile
}
## Decimo decil aparte para incluir el resto de observaciones
ing_med_de[n_percentil] <- sorted_ing %>% filter(cumulative_obs >= lim_inf) %>%
summarise(mean = sum(ing_cor*factor)/tam_decile) %>%
pull()
return(ing_med_de)
}
n_percentil <- 10
df <- tibble("2012" = ingreso_medio_percentil(concentradohogar2012, n_percentil),
"2014" = ingreso_medio_percentil(concentradohogar2014, n_percentil),
"2016" = ingreso_medio_percentil(concentradohogar2016, n_percentil),
"2018" = ingreso_medio_percentil(concentradohogar2018, n_percentil),
index = 1:n_percentil)
df <- df %>% gather(key = "Año", value = "Ingreso_medio", -index)
ingreso_medio_decil <- ggplot(df , aes(x = index, y = Ingreso_medio))+
geom_bar(stat="identity")+
facet_wrap(~ Año) +
scale_y_continuous(labels = scales::comma, breaks = seq(0,175000,25000))+
scale_x_continuous(breaks = 1:10, labels = scales::ordinal)+
xlab("Decil") + ylab("Ingreso medio") +
labs(title="Ingreso medio por deciles") +
theme_bw()
ingreso_medio_decilQuizá nos sea más fácil observar el comportamiento por decil si lo graficamos en términos de los porcentajes que acumula cada decil. Además esto nos dará una idea muy clara de la desigualdad presente entre el decil más rico y el decil más pobre.
ingreso_acumulado_percentil <- function(datos, n_percentil){
# Retorna el porcentaje del ingreso total que acumula cada percentil
# datos son las base de datos del censo, n_percentil es el numero de divisiones (e.g.: 10 retorna deciles, 100 retorna percentiles)
## Datos sorteados por ingresos trimestales corrientes más indice de los hogares representados observados
sorted_ing <- datos %>% select(ing_cor, factor) %>%
arrange(ing_cor) %>%
mutate(cumulative_obs = cumsum(factor))
## Total de hogares representados (ver definicion de factor en documentacion de la ENIGH)
total_hogares <- sum(datos$factor)
## Muestras por por percentil truncado para no incluir decimales
tam_decile <- trunc(total_hogares/n_percentil)
## Primeros n-1 deciles
lim_inf <- 1
lim_sup <- tam_decile
ing_total <- c()
for (i in 1:(n_percentil-1)) {
ing_total[i] <- sorted_ing %>% filter(cumulative_obs >= lim_inf, cumulative_obs <= lim_sup) %>%
summarise(mean = sum(ing_cor*factor)) %>%
pull()
## Limites del decil
lim_inf <- lim_inf + tam_decile
lim_sup <- lim_sup + tam_decile
}
## Decimo decil aparte para incluir el resto de observaciones
ing_total[n_percentil] <- sorted_ing %>% filter(cumulative_obs >= lim_inf) %>%
summarise(mean = sum(ing_cor*factor)) %>%
pull()
ing_total <- ing_total/sum(datos$ing_cor*datos$factor)
return(ing_total)
}n_percentil <- 10
df <- tibble("2012" = ingreso_acumulado_percentil(concentradohogar2012, n_percentil),
"2014" = ingreso_acumulado_percentil(concentradohogar2014, n_percentil),
"2016" = ingreso_acumulado_percentil(concentradohogar2016, n_percentil),
"2018" = ingreso_acumulado_percentil(concentradohogar2018, n_percentil),
index = fct_rev(as.factor(1:n_percentil)))
df <- df %>% gather(key = "Año", value = "Ingreso", -index)
ingreso_acu_decil <- ggplot(df, aes(fill = index, x = Año, y = Ingreso)) +
geom_bar(position="stack", stat="identity") +
ylab("Porcentaje de Ingreso acumulado") + xlab("Año") +
scale_fill_viridis(discrete = T, name = "Decil") +
scale_y_continuous(labels = scales::percent) +
theme_bw() +
labs(title= "Ingreso acumulado por Decil")+
theme(legend.background = element_rect(fill = "transparent", colour = "transparent"))
ingreso_acu_decilYa que esta ultima gráfica tampoco nos arroja información contundente sobre la consistencia de los datos a traves de los años, exploraremos qué ocurre con la acaparación de los ingresos en los ultimos cuantiles de los datos.
n_percentil <- 1000
df <- tibble("2012" = ingreso_acumulado_percentil(concentradohogar2012, n_percentil),
"2014" = ingreso_acumulado_percentil(concentradohogar2014, n_percentil),
"2016" = ingreso_acumulado_percentil(concentradohogar2016, n_percentil),
"2018" = ingreso_acumulado_percentil(concentradohogar2018, n_percentil),
index = 1:n_percentil)
df <- df %>% gather(key = "Año", value = "Ingreso", -index)
ingreso_acu_per <- ggplot(df, aes(x = index, y = Ingreso)) +
geom_line(aes(color = Año), size = 1) +
scale_x_continuous(labels = scales::ordinal) +
scale_y_continuous(labels = scales::percent) +
coord_cartesian(xlim=c(990,1000)) +
ylab("Porcentaje de Ingreso acumulado") + xlab("Percentil (0.1%)") +
labs(title="Ingreso acumulado por cada 0.1% de la poblacion") +
theme_bw() +
theme(legend.position = c(.1,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))
ingreso_acu_perOtra medida que nos puede arrojar información sobre posibles incosistencias es saber si hay gente que reportó que gasta más de lo que dijo que ingresaba. De la misma manera, lo haremos para los 4 censos bajo estudio para saber si podemos confiar en la métrica dada.
gastan_de_mas <- function(datos){
total_hogares <- sum(datos$factor)
datos %>%
transmute(mas_gasto = (gasto_mon > ing_cor)*factor) %>%
summarise(proporcion = sum(mas_gasto)/total_hogares) %>%
pull()
}
# Porcentaje de personas que gastan de más
over_spenders <- c(gastan_de_mas(concentradohogar2012),
gastan_de_mas(concentradohogar2014),
gastan_de_mas(concentradohogar2016),
gastan_de_mas(concentradohogar2018))
over_sp <- ggplot()+
geom_bar(aes(x=c("2012", "2014", "2016", "2018"), y=over_spenders), stat = "identity") +
scale_y_continuous(labels = scales::percent) +
xlab("Año") + ylab("")+
labs(title="Porcentaje personas que gastan más de lo que ingresan")+
theme_bw()
over_sp 20% es un porcentaje para nada despreciable de gente que al parecer tiene gastos corrientes mayores que sus ingresos. Lo que es más, está métrica es muy consistente en los 4 años por lo que es justo asumir que estamos frente a un claro fenomeno y no simple error humano.
Aquí veremos una forma de aporvechar los datos que nos brinda el SAT para mejorar nuestro entendimiento de los ingresos de las personas. En específico nos interesa saber cómo se comportan los ingresos más extremos que reportan frente al SAT. Por supuesto, a la hora de modelar habrá que hacer las suposiciones adecuadas para poder compararlo con aquel ingreso corriente que se reporta en la ENIGH. De esto se habla en el resumen que acompaña a este código.
breaks <- c(1,2,3,4,5,6,8,10, 15,20,30,50,100,200,500,1000)*1000000
means <- c()
for (i in 1:(length(breaks)-1)) {
means[i]<- mean(sat2012$I_DEC_TIAONCT1_AA[sat2012$I_DEC_TIAONCT1_AA >= breaks[i]
&
sat2012$I_DEC_TIAONCT1_AA<breaks[i+1]],
na.rm = T)
}
top_sat <- ggplot()+
geom_line(aes(x=1:(length(breaks)-1), y=means), color = "light blue", size = 1) +
scale_x_continuous(breaks = 1:(length(breaks)-1), labels = c("1-2 M", "2-3 M", "3-4 M", "4-5 M", "5-6 M", "6-8 M", "8-10 M", "10-15 M", "15-20 M", "20-30 M", "30-50 M", "50-100 M", "100-200 M", "200-500 M", "500-1000 M")) +
scale_y_continuous(labels = scales::comma, breaks=seq(0,600000000,100000000)) +
labs(title = "Ingreso medio segun ingresos reportados al SAT") +
xlab("Ingreso medio") + ylab("Rango de ingreso")+
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
top_satEn seguida les presento una forma para estimar parámetros con y sin restricciones. La especificación de las restricciones es muy importante para lograr una solución al problema de optimización. En este caso, solo busco presentarles el flujo de trabajo y no parámetros 100% reales a los presentados en los artículos de Bustos.
mean_log <- function(x) {
exp(x[1]+(x[2]^2)/2)
}
pll_log <- function(x){
-sum(dlnorm(concentradohogar2012$ing_cor+.01, x[1], x[2], log = TRUE)*concentradohogar2012$factor)
}
par_log_NR<- optim(par = c(10, 1), pll_log, control=list(maxit=1000))$par
par_log_R <- solnp(pars = c(11, 2), # Valores inciales
pll_log, # función a optimizar
eqfun=mean_log, # función restricción
eqB=92733.62, # Valor a igualar
LB=c(0,0), # Soporte mínimo para los parametros
UB=c(20,5))$pars # Soporte máximo para los parametros##
## Iter: 1 fn: 373023071.8492 Pars: 10.53334 1.83813
## Iter: 2 fn: 369313786.9753 Pars: 10.47885 1.57183
## Iter: 3 fn: 367073387.1403 Pars: 10.55500 1.36961
## Iter: 4 fn: 366549763.1526 Pars: 10.60202 1.29591
## Iter: 5 fn: 366512324.9051 Pars: 10.61045 1.28615
## Iter: 6 fn: 366511843.7034 Pars: 10.61068 1.28593
## Iter: 7 fn: 366511843.4866 Pars: 10.61068 1.28593
## Iter: 8 fn: 366511843.4875 Pars: 10.61068 1.28593
## solnp--> Completed in 8 iterations
## [1] 37999.64
## [1] 37254.33
## [1] 92733.62
Ahora observamos cómo se estima para otra distribución, en este caso la Beta generalizada de tipo II.
mean_gb2 <- function(x) {
moment.gb2(1, x[1], x[2], x[3], x[4])
}
pll_gb2 <- function(x){
-loglp.gb2(concentradohogar2012$ing_cor+.01, x[1], x[2],x[3], x[4], w=concentradohogar2012$factor)
}
par_gb2_NR<- optim(par = c(1.307071, 20594.790235, 2.341514, 1.843004), pll_gb2, method = "BFGS",control=list(maxit=1000))$par
par_gb2_R <- solnp(pars = c(1.5, 30000, 1, 1), # Valores inciales
pll_gb2, # función a optimizar
eqfun=mean_gb2, # función restricción
eqB=92733.62, # Valor a igualar
LB=c(0,0,0,0), # Soporte mínimo para los parametros
UB=c(5,100000,5,5))$pars # Soporte máximo para los parametros##
## Iter: 1 fn: 11.5445 Pars: 1.47087 30642.59739 1.32305 1.01227
## Iter: 2 fn: 11.5258 Pars: 1.48804 30640.46754 1.04509 0.94058
## Iter: 3 fn: 11.5258 Pars: 1.48801 30640.47383 1.04545 0.94067
## Iter: 4 fn: 11.5258 Pars: 1.48801 30640.47383 1.04545 0.94067
## solnp--> Completed in 4 iterations
## [1] 37999.64
## [1] 38249.42
## [1] 92733.62
Ahora nos gustaría saber qué tanta corrección hubo con respecto de la distribución original que nos arrojan los datos. Es decir, cómo afecta a las distribuciones estimadas el uso de una restricción en su media. De nuevo, una plática más extensa se encuentra en el resumen.
distribucion_empirica <- function(x, datos, factor){
# x es el valor donde se quiere evaluar la función de distribución empírica
# datos es el conjunto de datos en formato de vector que se quiere evaluar
df <- data.frame(datos = datos, factor = factor, cumu = cumsum(factor))
F_x <- c()
for (i in 1:length(x)) {
F_x[i] <- sum(df[df$datos<x[i], "factor"])/df[length(datos), "cumu"]
}
return(F_x)
}
distri_empi_restri <- ggplot() +
geom_function(aes(colour = "Empirica"),
fun = distribucion_empirica,
args = list(datos = concentradohogar2012$ing_cor,
factor = concentradohogar2012$factor),
size = 1) +
geom_function(aes(colour = "GB2"),
fun = pgb2,
args = list(shape1 = par_gb2_R[1],
scale = par_gb2_R[2],
shape2 = par_gb2_R[3],
shape3 = par_gb2_R[4]),
size = 1) +
geom_function(aes(colour = "Log-Normal"),
fun = plnorm,
args = list(meanlog = par_log_R[1],
sdlog = par_log_R[2]),
size = 1) +
scale_x_continuous(trans = 'log', limits = c(600, 4000000), breaks = 10^(seq(2,6,1)), labels = scales::comma) +
scale_y_continuous(breaks = seq(0,1,.1)) +
scale_colour_manual(name="Distribución",values=c("Empirica"="#3591d1", "GB2"="#f04546", "Log-Normal"="#f045d6")) +
ylab("Probabilidad") + xlab("Ingreso (Escala log)") +
labs(title = "Distribución empírica contra estimadas con restricción") +
theme_bw() +
theme(legend.position = c(.1,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))
distri_empi_restriSiguiendo el mimso proceso pero ahora presento las distribuciones con y sin restricciones. Vemos que las correciones son más o menos consistentes entre los modelos. Esto son buenas noticias pues podemos hacer algunas inferencias sin preocuparnos de casarnos con una distribución en específico.
p1 <- ggplot() +
geom_function(aes(colour = "Log-Normal Restringida"),
fun = plnorm,
args = list(meanlog = par_log_R[1],
sdlog = par_log_R[2]),
size = 1) +
geom_function(aes(colour = "Log-Normal No Restringida"),
fun = plnorm,
args = list(meanlog = par_log_NR[1],
sdlog = par_log_NR[2]),
size = 1) +
scale_x_continuous(trans = 'log', limits = c(600, 4000000), breaks = 10^(seq(2,6,1)), labels = scales::comma) +
scale_y_continuous(breaks = seq(0,1,.1)) +
scale_colour_manual(name="Distribución",values=c("Log-Normal Restringida"="#f04546", "Log-Normal No Restringida"="#3591d1")) +
ylab("Probabilidad") + xlab("Ingreso (Escala log)") +
theme_bw() +
theme(legend.position = c(.25,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))
p2 <- ggplot() +
geom_function(aes(colour = "GB2 Restringida"),
fun = pgb2,
args = list(shape1 = par_gb2_R[1],
scale = par_gb2_R[2],
shape2 = par_gb2_R[3],
shape3 = par_gb2_R[4]),
size = 1) +
geom_function(aes(colour = "GB2 No Restringida"),
fun = pgb2,
args = list(shape1 = par_gb2_NR[1],
scale = par_gb2_NR[2],
shape2 = par_gb2_NR[3],
shape3 = par_gb2_NR[4]),
size = 1) +
scale_x_continuous(trans = 'log', limits = c(600, 4000000), breaks = 10^(seq(2,6,1)), labels = scales::comma) +
scale_y_continuous(breaks = seq(0,1,.1)) +
scale_colour_manual(name="Distribución",values=c("GB2 Restringida"="#f04546", "GB2 No Restringida"="#3591d1")) +
ylab("Probabilidad") + xlab("Ingreso (Escala log)") +
theme_bw() +
theme(legend.position = c(.2,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))
grid.arrange(p1,p2,ncol=2)Por último, revisemos una muestra rápida de cómo se ve la desigualdad estimada cuando asumimos que los modelos antes descritos son los correctos frente a la desigualdad que existe entre la muestra del ENIGH.
## Calculamos el acumulado por decil mediante simulación
x <- rlnorm(100000, par_log_R[1], par_log_R[2])
deciles <- quantile(x, probs = seq(0,1,.1))
total<-sum(x)
medd <- c()
for (i in 1:10) {
medd[i] <- sum(x[x > deciles[i] & x <= deciles[i+1]])/total
}
## Construcción de la tabla para alimentar el ggplot
n_percentil <- 10
df <- tibble("ENIGH 2012" = ingreso_acumulado_percentil(concentradohogar2012, n_percentil),
"Log-Normal" = medd,
index = fct_rev(as.factor(1:n_percentil)))
df <- df %>% gather(key = "Fuente", value = "Ingreso", -index)
## Construcciónd de la gráfica
ingreso_acu_decil_ENIGH_Log <- ggplot(df, aes(fill = index, x = Fuente, y = Ingreso)) +
geom_bar(position="stack", stat="identity", width=0.7) +
ylab("Porcentaje de Ingreso acumulado") + xlab("Modelo") +
scale_fill_viridis(discrete = T, name = "Decil") +
scale_y_continuous(labels = scales::percent) +
theme_bw() +
theme(legend.background = element_rect(fill = "transparent", colour = "transparent"))
ingreso_acu_decil_ENIGH_LogPodemos observar que la desigualdad se acrecenta considerablemente. Esto por que el decil más rico acapara un porcentaje aún mayor de el ingreso corriente.
Se generan los individuos y se especifica el total del dinero en la simulación
Se hace un dataframe para guardar la información de las simulaciones
Se realiza la simulación
for (i in 1:(4*10000)){
# elegir los dos agentes de la transacci?n
primer_agente = round(runif(1, 1, 500))
segundo_agente = round(runif(1, 1, 500))
while (primer_agente == segundo_agente){
segundo_agente = round(runif(1, 1, 500))
}
# determinar el ganador
ganador = rbinom(1, 1, 0.5)
if (ganador == 1){
indice_ganador = segundo_agente
indice_perdedor = primer_agente
} else {
indice_ganador = primer_agente
indice_perdedor = segundo_agente
}
# calcular el dinero a ser transferido usando la regla tres
transaccion = runif(1, 0, 1) * din_total/N
# checar si el perdedor tiene suficiente dinero para que se
# realice la transacci?n
if (dataind[indice_perdedor, "dinero"] < transaccion){
next
}
# llevar a cabo la transacci?n
dataind[indice_ganador, "dinero"] = dataind[indice_ganador, "dinero"] + transaccion
dataind[indice_perdedor, "dinero"] = dataind[indice_perdedor, "dinero"] - transaccion
}Una vez que se ha realizado la simulación, se grafica un histograma para obtener la densidad de la distribución del ingreso para los agentes y sobre este histograma, se grafica la función de densidad propuesta para modelar la distribución de este ingreso
hist(dataind[,"dinero"],freq = FALSE)
lines(seq(1:7000), (densidad = (1/1000) * exp(-seq(1:7000) / 1000)))